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Abstract 

A novel development is given of the theory of Gaussian quadrature, not relying 
on the theory of orthogonal polynomials. A method is given for computing the 
nodes and weights that is manifestly independent of choice of basis in the space 
of polynomials. This method can be extended to compute nodes and weights for 
Gaussian quadrature on the unit circle and Gauss type quadrature rules with 
some fixed nodes. 



The aim of this letter is to show how the theory of Gaussian quadrature can be developed 
using elementary linear algebra, without relying on the theory of orthogonal polynomials. 
This is certainly useful from a pedagogical point of view, showing how basic results in linear 
algebra can have powerful applications and allowing Gaussian quadrature to be taught in 
introductory courses before the rather heavier topic of orthogonal polynomials. But it is 
also important from a more fundamental point of view. Orthogonal polynomials are an 
invaluable technical tool, but they are, after all, just a choice of basis in the vector space 
of polynomials. In this paper we give a method for computing nodes and weights that is 
independent of the choice of basis. 

A number of the formulae in this letter can also be found in the papers of Sack and 
Donovan [7] and Gautschi [5] on the "modified moment method" for computing nodes and 
weights. The emphasis, however, is somewhat different: these papers focus on efficient, sta- 
ble computation, while the focus here is more on the conceptual issue of basis indepedence. 
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It also turns out that occasional benefit can be obtained by a choice of basis not covered in 
[7] and [5]. Our basis independent method for computing nodes and weights is described in 
theorem 3. It allows immediate generalization to the cases of Gaussian quadrature on the 
unit circle and Gauss type quadrature with some fixed nodes, as we shall show. Our view 
of quadrature formulae as being related to blinear forms allows a generalization to multiple 
dimensions, though this is discussed in a separate paper [3]. 

There is, of course, already a completely elementary approach to Gaussian quadrature 
that avoids mention of orthogonal polynomials. We seek Xi,Wi, i = 1, . . . , n + 1, such that 
the quadrature formula 

rb n+l 

/ w(x)f(x)dx « W if( X i) , Wi^O , (1) 
Ja i=l 

is exact when f(x) is any polynomial of degree not more than 2n + 1. Here w(x) is a 
known non-negative weight function on the interval [a, b]. Substituting in turn f(x) = 
1, x, x 2 , . . . , x 2n+l we see that the nodes and weights must satisfy the equations 

n+l 

J2 w i x i= m a, a = 0, 1, . . . ,2n + 1 , (2) 

i=l 



where m a is the ath moment 

m 



rb 

a = w(x)x a dx . (3) 

J a 



In principle we should be able to obtain the nodes and weights by solving the system (direqs) 
of nonlinear equations. The drawback of this approach is that there is no evident way to 
solve (direqs) except in very simple cases; it is far from clear that there even exists a solution. 
Also, dependence on choice of basis is hardly resolved in this approach. 

A linear algebra approach to Gaussian quadrature. Let V n denote the vector space 
of polynomials of degree up to n (so dimV n = n + l). Let w(x) be a given non-negative 
weight function on the interval [a, b], and write 

rb 

< f\g >= / w(x)f(x)g(x) dx (4) 

J a 

for any polynomials /, g. For any n, < -\- > determines an inner product on the space V n . 
Define the symmetric bilinear form X(-, •) on V n by 

rb 

X(f,9)= w(x)xf(x)g(x) dx . (5) 

J a 

For any symmetric bilinear form on an inner product space there exists an orthonormal basis 
in which the form is diagonal. In other words, we can find eigenvalues x±, . . . , x n+ i G R and 
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eigenvectors Ui(x), . . . , u n+1 {x) G V n such that 

<Ui\uj > = 5ij , (6) 
X(ui,Uj) = Xi5ij . (7) 

Theorem 1. Let x\,..., x n+ i be the eigenvalues of the bilinear form X(-,-) with correspond- 
ing orthonormal eigenvectors Ui(x), . . . , u n+ i(x). Then for any polynomial p(x) G Vi n +\ 

rb n + 1 

/ w(x)p(x)dx =^2 < l\ui> 2 p(xi) , (8) 
Ja i=i 

i.e. the quadrature rule 

rb n + 1 

/ w(x)f(x)dxttJ2wif(xi), (9) 

where Wi =< \\ui > 2 , is exact whenever f(x) G Vi n +\- 
The proof proceeds via a lemma. 

The S Lemma. For any polynomial p G V n +i, we have 

< p\ui >=< l\ui > p(x { ) . (10) 

In other words, up to a normalization factor, the projection of p(x) on Ui(x) is found by 
simply evaluating p(x) at x — £3. 

Proof of The 5 Lemma. We prove the result for p(x) = 1, x, . . . , x n+1 , the general result 
follows by linearity. For p(x) = 1 the result is trivial. For j — 1, . . . , n + 1 we have 

rb , 

< X^\Ui >= / -U^X^Ii^x) <i£ = X(x- 7 ^ 1 ,Mj) = Xj < X J_1 |'Uj > . (11) 
J a 

Iterating this procedure j times we obtain < x^\ui >= x\ < l\ui >, as required. • 

Proof of Theorem 1. We prove the result for p(x) = l,x, . . . ,x 2n+1 , the general result 
follows by linearity. For p(x) = 1 we have 

/b n+1 n+1 

w(x) dx =< 1 1 1 >=J2< l \ Ui >< ^l 1 >= H < l \ Ui >2 ( 12 ) 
i=i i=i 

as required. For p(x) = x\ j — 1, . . . , 2n + 1, we can write p(x) = x J ' 1+J ' 2+1 , where ji, j'2 G 
{1,2,..., n}, and thus 

w(x)x 1 dx = X(x n ,x n ) 

n+1 n+1 

= J2J2 < xJ1 \ u h >< xn \ U ii > X (u h ,u i2 ) 

h=l 12=1 
n+1 

= ^<ii^> 2 xi, (13) 
i=i 
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where in the last step we have used the 5 Lemma. • 

Theorem 1 proves the existence of a Gaussian quadrature formula for any interval [a, b] 
and non-negative weight function w(x). The proof is based just on diagonalization of the 
bilinear form X(-, •) with no prior knowledge needed of orthogonal polynomials. The next 
theorem shows that some of the main results on Gaussian quadrature are easily obtained 
within our approach. 

Theorem 2. 

1) The Gaussian quadrature formula of theorem 1 is unique. 

2) The nodes of the Gaussian quadrature formula lie in the interval [a, b]. 

3) The nodes of the Gaussian quadrature formula are roots of any nontrivial degree n + 1 
polynomial that is orthogonal to all polynomials of degree at most n. 

To motivate the first part of the proof we make the observation that by setting p(x) = Uj(x) 
in the 5 Lemma we obtain 



implying the interesting result that the functions Ui(x) are actually proportional to the 
Lagrange cardinal functions for the points x±, . . . ,x n+ i. Note that by setting p(x) = Ui(x) 
in (dl) we see that < l\ui > cannot be zero. Note also that the last equation implies the 
Xi are distinct. We mention in passing that the fact that Ui(xj) is proportional to 5^ is one 
reason we gave the 5 Lemma its name. 

Proof. 1) Suppose X i: W i: i = 1, . . . , n + 1 are nodes and weights for a Gaussian quadrature 
formula, i.e. for any polynomial p(x) G V2n+i 



with the Xi distinct and Wj > 0. Define the nth degree polynomials Ui(x), . . . , U n+ i(x) by 
requiring 




(15) 





n + 1 . 



(16) 



We then have 



J a 



b 



n+1 



< Ui\Uj > 



w(x 



Wi^Ujix) dx = Y J W k Ui{X k )Uj{X k ) = 5, 



(17) 



k=i 



and 




n+1 



w(x)xU i (x)U j (x) dx=Y, W k X k Ui{X k )Uj{X k ) = X&j . 



(18) 



k=i 
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Thus the X { are necesarilly the eigenvalues of the bilinear form X, and the Ui are the 
orthonormal eigenvectors. Furthermore we have 

rb n + l , 

< l\Ut >= / w(x)Ui(x) dx=J2 W k Ui(X k ) = ^Wi , (19) 
Ja k=i 

so Wi =< l\Ui > 2 , as before. 

2) We have 

rb rb rb 

/ w(x)au 2 (x) dx < w(x)xu 2 (x) dx < w{x)bu 2 (x) dx , 

J a J a J a 

so a < X(ui, Ui) < b and therefore a < Xi < b. 

3) Note that the 5 Lemma holds for all p(x) £ V n +i- If p(x) is a polynomial of degree n + 1 
that is orthogonal to all polynomials of degree at most n, then it is orthogonal to all the 
Ui(x), and thus, by the 8 Lemma, p(xi) = for all i. • 

Calculation of the nodes and weights. Having presented our nonstandard development 
of the theory of Gaussian quadrature, we turn our attention to the question of computing 
nodes and weights using an algorithm that allows for an arbitrary choice of basis in V n . To 
find the nodes and weights we clearly need information on the inner product < -|- > and 
the blinear form X(-, •) on V n . So assume the functions qi(x), q 2 (x), . . . , q n+ i(x) are a basis 
of V n and that we are given the (n + 1) x (n + 1) matrices B, A defined by 

rb 

Bij = < qMj >= / w{x)q i {x)q j {x)dx , i, j = 1, . . . ,n + 1 , (20) 

J a 

rb 

= X(q i ,q j )= w(x)xq i (x)q j (x)dx , i, j = 1, . . . , n + 1 . (21) 

J a 

We are about to see that knowledge of the matrices B, A and just one of the functions qi(x) 
is sufficient to determine the nodes Xi and weights Wi. As a prelude, note that there exists 
a diagonal matrix D and a matrix V satisfying 

V T BV = I , AV = BVD . (22) 

This is nothing but the statement that X(-,-) has orthonormal eigenvectors, written in 
the basis In Matlab Release 14 the matrices D, V are computed from A, B by the 

command [V D]=eig(A,B, 'chol' ). 

Theorem 3. If D, V are matrices satisfying (VD), with D diagonal, then the nodes and 
weights of the Gaussian quadrature formula (that is exact for polynomials in T^n+i) are 
determined by 

x t = Da , Wi = (^r¥) ■ (23) 
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(The second relation holds for arbitrary j). 

Proof. Since all the integrals needed to compute B and A are integrals of polynomials of 
degree up to 2n + 1, we have 

n+1 n+1 

B ij = J2 VJ k ( li( X k)(!j(Xk) = ^2QikW k kQjk , (24) 
k=l k=l 
n+1 n+1 

A v = J2 w k x kQi(xk)qj(x k ) = ^2QikWkkX k kQjk , (25) 

k=l k=l 

where we have introduced diagonal matrices X, W with the nodes X{ and the weights Wi 
along their diagonals, and the matrix Q with entries Qik = qi(xk)- More concisely, we have 

B = QWQ T , A = QWXQ T . (26) 

Inserting these into the second relation in (VD) we obtain 

QWXQ T V = QWQ T VD , (27) 

and therefore 

X = (Q T V)D(Q T Vy 1 . (28) 

Since diagonal matrices are conjugate if and only if they are equal, up to reordering of 
the diagonal entries, we deduce the first statement in the theorem. Furthermore, after 
reordering the nodes if necessary, we see that the matrix Q T V commutes with X. But X 
is a diagonal matrix with distinct entries on the diagonal, so Q T V must also be diagonal. 
Write VWQ T V = Y. Then Y T Y = V T QWQ T V = I (by (rheqs) and the first relation in 
(VD)). So the diagonal entries of Y are all plus or minus 1. Writing y/W Q T = YV^ 1 we 
obtain WiQ 2 ^ = ((V -1 )^) 2 for all i,j, giving the second statement in the theorem. • 

Notes: I. Under a change of the basis {qi} we will have Q — > MQ, A — > MAM T , 
B — > MBM T , V — > M~ T V, where M is a suitable nonsingular matrix. The nodes and the 
weights, however, remain unchanged. 

2. Once the nodes Xj are known, then the weights can be obtained from the first equation 
in (rheqs). This, however, requires full knowlegde of the basis g«, whereas the formula for 
the weights in the theorem requires knowledge of just a single basis element. 

3. The simplest cases of theorem 3 are well-known. If we choose qi(x) to be the orthormal 
basis of V n formed by Gram-Schmidt orthonormalization of the basis l,x, ...,x n , then 
B = I and A is tridiagonal. The method for finding nodes and weights reduces to the 
classic and widely used method of [6] (see [1] for an exposition). If we choose (&(x) = x % ~ x 
then the matrices A and B are Hankel matrices (Aij and B^ both depend only on i+j), and 
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the computation of D and V is numerically unstable (at least for all but the smallest values 
of n). For more general bases the method is essentially equivalent to that of [7] and [5], but 
the formulation here makes the issue of basis independence much clearer. The papers [7] 
and [5] focus on bases in which the g«(x) satisfy a recursion of the form 



xqi(x) = aiq i+1 (x) + hq^x) + dq^x) 



(29) 



This makes most of the entries of A expressible in terms of the entries of B and the coeffi- 
cients aj, 6j, Cj. 

Example: For the weight function w(x) = (1 + x)" 1 on [0, 1], the construction of orthog- 
onal polynomials by Gram-Schmidt orthonormalization of the standard basis 1, x, . . . , x n is 
awkward. It is much more convenient to work in a basis in which most of the elements have 
a factor (1 + x). As a first attempt, take 



qi(x) 



(1 + x)x l 1 i — 1, . . . , n 
1 i = n+ 1 



(30) 



It is straightforward to compute the necessary integrals to obtain 



Bij = < 



i 



i+j-i 



i 



i 

i+j 



1 

3 

In 2 



i,j = l,...,n 
i — 1, . . . ,n , j = + l 
j — 1, . . . ,n ,i — n + 1 

i — j — n + 1 



(31) 



r — + - J — 

1 



i+1 
1 

i+i 
1 - In 2 



1,.. 



: 1, . . .,n 
n , j — n + 1 



(32) 



j = l,...,n,z = n+ l 
i = j = n + 1 

For fairly small n, it is straightforward to find D, V and the nodes and weights using Matlab. 
As n gets larger and the conditioning of the Hankel matrices B and A becomes a concern, 
it is preferable to use the basis 



q%{x) 



[l + x)pi(x) 
1 



i— 1, . . . , n 
i = n + 1 



(33) 



where the {pi(x)} are the standard orthonormal polynomials on [0,1] (orthonormal with 
respect to the weight function 1). Then B and A are tridiagonal and pentadiagonal matrices 
respectively. 

The point in this example is that we can exploit the freedom of basis to simplify evaluation 
of the necesary integrals. The cost is a little more linear algebra, but this is hardly a challenge 
by modern standards. 
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r-2n . . "+ 1 

/ / (e id ) w(6)d9 » ]T ^/(^ , w r ^ , (34) 



Generalization 1. Gaussian Quadrature on the Unit Circle [2]. Let w(0) be a 
2n— periodic non-negative weight function. Assume there exists a quadrature formula 

1 r 2ir , . "+ 1 

r=l 

with distinct nodes z r , which is exact for the In + 2 functions /(z) = £~ ra , . . . , z n , z ra+1 and 
all their linear combinations. Let q r (z), r = 1, . . . , n + 1 be a basis for the vector space of 
(complex) polynomials of degree at most n in z. Define the (n + 1) x (n + 1) matrices S, A 
by 

1 r2ir 

B rs = — / w(0)g r (e ifl ) gs (e lfl )^ , r, s = 1, . . . , n + 1 , (35) 

27T JO 

A rs = — / ^% r (e^yV(e*V0 , r,5 = l,...,n + l. (36) 

27T JO 

Note that 5 is Hermitian and positive definite, but A need not be Hermitian. We assume 
the pair A, B is diagonalizable, i.e. that there exists a diagonal matrix D and a matrix V 
such that 

AV = BVD . (37) 

The matrices D, V can be obtained in Matlab by the command [V D]=eig(A,B). We then 
have the following analog of theorem 3: 

Theorem 3'. The nodes and weights of the quadrature formula (quad3) (that is exact for 
f(z) = z~ n , . . . , z n , z n+1 ) are determined by 

z r = D rr , w r = ^=^=- . (38) 

q s (z r )q s (i) 

(The second relation holds for arbitrary s). 
Proof: Proceeding as in the real case we obtain 

B = QWQ , A = QWXQ , (39) 



where Q, Q are matrices with entries Q rs = q r (z s ) and Q rs = q s (=), respectively, and X, W 
are diagonal matrices with diagonal entries equal to the nodes and the weights, respectively. 
Using these relations in (VD2) we rapidly obtain 

X = (QV)D(QV)- 1 , (40) 

giving us the first result, that X = D, and also that the matrix Y = QV is diagonal. 
From Q = YV^ 1 we deduce that 

Qrs l^. r (V^ ) rs , (41) 



8 



for each r, s. Multiplying the first equation in (bqw) on the right by V, gives BY = QWY, 
implying 

(BV) sr = Y rr w r Q sr , (42) 

for each r, s. Finally we divide (tyl) by (ty2) to obtain the second result in the theorem. • 

Notes: 1. The result in the theorem shows how to compute the weights given A, B and 
one of the q r . It may be easier, if all of the q r are known, to directly use the first formula 
in (bqw) to find the weights. 

2. For each k we have Av k = z k Bv k , where v fc denotes the vector that is the A;th column of 
V. Thus v* k (A - z k B)v k = 0, or 

r2ix 



2tt Jo 



(e i9 - z k )\a(9)\ 2 w(9)d9 = , (43) 

Z7I JO 

where a(9) = I] r (v fc ) r g r (e* 6 '). Since 

/ e^ l9 \a(9)\ 2 w(9)d9 < \a(9)\ 2 w(9)d9 , (44) 
Jo Jo 



it follows that \z k \ < 1, i.e. the nodes are inside the unit circle. 

Example: Taking w{9) = sin 2 9 and using the standard basis {q r (z)} = {1, z, . . . , z n } we 
have 

1 r 2n , 1 

A rs = — e^ 9 sin 2 9d9 = - (25 r - s - 5 r _ s+2 - 5 r _ s _ 2 ) , (45) 

2 71 jo 4 

B rs = ±- r e i( - r - s+i y> sin 2 9d9 = - (25 r _ s+1 - 5 r _ s+3 - 6^) , (46) 
1-k Jo 4 

where 5 r is 1 if r = and otherwise, and we have used ^ / 27r e* r£ '(i6' = 5 r . With n = 7 
and these simple choices of ^4 and B, and with a little help from Matlab, theorem 3' rapidly 
reproduces the nodes and weights given in the third part of table 1 in [2]. 

For a further reference on the use of Gaussian quadrature for complex integrals see [8]. 

Generalization 2. Gauss-type rules with some fixed nodes. We seek nodes and 
weights Xi,Wi, i = 1, . . . , n + 1, and weights v a , a = 1, . . . , m, such that the quadrature 
formula 

/fe m n+l 

f(x)w(x)dx W VafiVa) + W if( X i) j V a ^ , (47) 

a=l i=l 

is exact when f(x) is any polynomial of degree not more than 2n + m + 1. Here the y a , 
a = 1, . . . , m, are given nodes in [a, b], and we assume the x^ and y a are all distinct. There 
is no particular reason to assume w(x) is non-negative. Let qi(x), i = 1, . . . ,n + 1 be an 
arbitrary basis for the vector space of polynomials of degree at most n. Let the (n+l) x (n+l) 
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matrices B, A be defined by 

Va) qi{x)qj{x)dx , i, j = 1, . . . ,n + 1 , (48) 

y a ) xqi(x)qj(x)dx , i, j = 1, . . . ,n + 1 . (49) 

Both B and A are symmetric, but in this generalization B need not be positive definite 
(though in the important cases m — 1, y\ — a or b, and m — 2, yx — a, y% — b, it is 
positive or negative definite). Assume the pair A, I? is diagonalizable, i.e. that we can find 
a diagonal matrix D and a matrix V satisfying 

AV = BVD . (50) 

We then have the following: 

Theorem 3". If the quadrature formula (quad2) is exact when f(x) is any polynomial of 
degree not more than In + m + 1, and the pair A, B is diagonalizable, then the nodes and 
weights are determined by 

x t = Da , wi ft { Xi - y a ) = (X^M^pii . (51) 

a=l Wi) 

(The second relation holds for arbitrary j). 

The proof of this is sufficiently similar to that of theorem 3' that we omit it. 

Once the nodes X{ and weights Wi have been determined using this theorem, the weights 
v a can be determined by solving the linear system arising from exactness of the quadrature 
formula for f(x) — 1, x, . . . , x m ~ x . 
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